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We consider asymmetric spin-1/2 two-leg ladders with non-equal antiferromagnetic (AF) couplings 
J|l and kJu along legs (k < 1) and ferromagnetic rung coupling, J±. This model is characterized 
by a gap A in the spectrum of spin excitations. We show that in the large J± limit this gap is 
equivalent to the Haldane gap for the AF spin-1 chain, irrespective of the asymmetry of the ladder. 
The behavior of the gap at small rung coupling falls in two different universality classes. The first 
class, which is best understood from the case of the conventional symmetric ladder at k = 1, admits 
a linear scaling for the spin gap A ~ J±. The second class appears for a strong asymmetry of the 
coupling along legs, kJ|| <C J± <C J\\ and is characterized by two energy scales: the exponentially 
small spin gap A ~ J± exp(— Jy / J±), and the bandwidth of the low-lying excitations induced 
by a Suhl-Nakamura indirect exchange ~ j]_/J\\. We report numerical results obtained by exact 
diagonalization, density matrix renormalization group and quantum Monte Carlo simulations for 
the spin gap and various spin correlation functions. Our data indicate that the behavior of the string 
order parameter, characterizing the hidden AF order in Haldane phase, is different in the limiting 
cases of weak and strong asymmetry. On the basis of the numerical data, we propose a low-energy 
theory of effective spin-1 variables, pertaining to large blocks on a decimated lattice. 



PACS numbers: 75.10.Pq, 71.10.Fd, 73.22.Gk 

I. INTRODUCTION 

Recent progress in nanotechnologies, molecular elec- 
tronics and quantum computing reinvigorated the inter- 
est to low-dimensional systems. Special attention has 
focused during recent years on quantum dots, arrays of 
coupled quantum dots, 1 quantum wires, spin chains or 
ladders. Another class of physical systems where low- 
dimensionality can be achieved is ultra-cold gases in op- 
tical lattices, 2 ' 3 which form good prototype systems for 
investigation of many strongly-correlated effects, such as 
metal-insulator transition, low-dimensional superconduc- 
tivity or formation of various density-wave states. The 
advantage of these systems is the high controllability of 
model parameters with external fields and preparation 
conditions. 

Many of these systems display low-energy feature that 
fall outside the standard behavior predicted by Landau's 
Fermi liquid or symmetry breaking theories. In par- 
ticular, the existence of a non-local (string) order pa- 
rameter is proven, both analytically 4 and numerically 5 , 
to be a characteristic feature of several classes of 
one-dimensional (ID) and quasi-lD systems. 6 ' 7 Such 
non-local order parameters are topologically protected 
against any local perturbations. The nature of these 
order parameters and the connection to topological 
invariants 6 is well understood 8 for spin chains with 
S > 1. For instance, the Haldane conjecture 9 proposed 
more than 20 years ago states that the properties of 
SU(2) symmetric antiferromagnetic (AF) spin-S 1 Heisen- 



berg chains differ for integer and half-integer spins. The 
excitations in the AF Heisenberg chains with half-integer 
spins are gapless 10 whereas in the integer spin case, a gap 
is present. The pure one-dimensional (ID) AF spin-1/2 
Heisenberg chain can be mapped onto a Luttingcr liquid 
which allows an exact bosonization treatment, resulting 
in a well understood gapless phase. 10 In contrast, for the 
AF spin-1 Heisenberg chain it is widely accepted that 
the excitations exhibit a gap, thanks to extensive numer- 
ical n ~ 16 and experimental 17,18 analysis. 

The present understanding of systems of coupled iden- 
tical 5 = 1/2 chains (spin ladders) is based on its similar- 
ities to larger-spin chains. This similarity allows the one- 
to-one translation of Haldane's conjecture, originally for- 
mulated for large-spins chain, onto spin-1/2 ladders with 
an odd or even number of identical legs. The usual as- 
sumption about the equivalence of the individual chains 
constituting the ladder, is referred below as a symmetric 
ladder situation. While the behavior of symmetric lad- 
ders has been thoroughly investigated, both theoretically 
and experimentally 6 ' 19 , the case of spin ladders with in- 
equivalent legs (asymmetric ladders) is less understood. 
The behaviour of the spin gap, in particular for the case 
of a single chain coupled to nearly free spins ("dangling 
spins") was recently discussed in 20-23 , but no firm conclu- 
sions about the gap scaling in the weak coupling regime 
were made there. 

In this work, we consider the Spiral Staircase Heisen- 
berg ladder (SSHL), consisting of of two unequal antifcr- 
romagnetically coupled spin-1/2 chains with a ferromag- 
netic (FM) rung coupling J±. Geometrically, this model 
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FIG. 1: (Color online) (a) Sketch of the Spiral Staircase 
Heisenberg Ladder, (b) View of the system from the top. 
For 8 — the model corresponds to the standard antifcr- 
romagnetically coupled spin- 1/2 Heisenberg ladder with FM 
rung coupling. The case 9 — n corresponds to the ID ST/ (2) 
symmetric single-pole ladder model. 



may be understood as a continuous twist deformation of 
an isotropic 2-leg ladder with interleg coupling Jy along 
leg 1 by an angle 9 (see Fig. 1). As a result of the de- 
formation, the coupling between neighboring sites along 
leg 2 is rescaled to the form Jy cos 2 (6/2) leading to the 
Hamiltonian (1) below. 

As is known from the seminal bosonization work of 
Shelton, Nersesyan and Tsvelik, 4 the FM coupling of AF 
spin chains induces a Haldane gap proportional to J± for 
weak interleg couplings. However, when the spin veloc- 
ity on one of the legs vanishes the bosonization method 
fails as seen from the fact that a simple formulation of 
the continuum limit, on which the bosonization approach 
relies, is inhibited. Alternative approaches, such as us- 
ing other analytical or numerical methods, are then de- 
manded. Beside the theoretical interest, an experimental 
motivation comes from the fact that the single-pole lad- 
der model at 9 = it can been used for modeling a stable 
organic biradical crystal PNNNO. 24 

In our previous work (Ref. 23), we used the quantum 
Monte Carlo simulations for investigation of asymmetric 
ladders. We demonstrated a non-zero value of string or- 
der parameter for the whole family of such ladders, which 
confirmed that the system is in a Haldane phase. We 
also presented numerical evidence for the smaller energy 
scale, J±/ J\\, associated with Suhl-Nakamura interaction 
(see below). Numerical results for the spin gap were also 
judged compatible as vanishing as J±/J\\, even though a 
faster decay could not be ruled out. These results were 
consistent with the flow equation calculations of Essler 
and coworkers 22 . 

In the present paper we show, that the spin gap 
changes its behavior at small J± <C Jy from A ~ J± in 
the symmetrical ladder case to A ~ J± exp(— J||/Jj_) in 
the single-pole ladder case. Alternatively, we may think 
about rungs with a fixed exchange coupling Jj_, while 



the leg exchange increases gradually. We observe that 
the small- Jy behavior A ~ Jy is the same for all asym- 
metric ladders, but that important differences appear in 
the limit J|| — > oo. The symmetrical ladder displays a 
saturation of the spin gap, A ~ Jj_, while the gap in 
the single-pole ladder reaches a maximum at Jn ~ Jj_ 
and exhibits exponential suppression beyond that scale, 
according to the above formula. 



The plan of the paper is as follows. We introduce the 
model and important definitions in Section II, where we 
also provide simple qualitative considerations. Particu- 
larly, we explain here the importance of Suhl-Nakamura 
(SN) indirect exchange between dangling spins in the 
single-pole situation. 



We propose analytical approaches to our model in Sec- 
tion III. In Section III A we develop a theory, incorpo- 
rating spin-wave analysis, SN interaction and Kadanoff 's 
decimation procedure, which satisfactorily describes the 
whole body of numerical data presented below. We stress 
that several quantities are extracted from the numerical 
data and plotted specifically for comparison with the pre- 
dictions of this effective theory. The unusual slow satu- 
ration of the spin gap value at large J± is discussed in 
Section IIIB. 



In Section IV we describe the results of our extensive 
numerical investigations of the problem. Section IV A 
discusses exact diagonalization (ED) results. The ap- 
pearance of a SN energy scale J]_/J\\ is shown here. Given 
the long-range character of the SN interaction which re- 
quire large systems to be appreciated, we resort to large 
scale quantum Monte Carlo (QMC) simulations in the 
next subsection. With this approach, we investigate the 
form of the spin correlation functions in Section IV B, 
compute the spin gap in Section IV C as well as the string 
order parameter in Section IV D, characteristic of the 
Haldane phase. The most challenging case of the single- 
pole ladder gives rise to the largest uncertainties for the 
spin gap, and we focus on this system with the use of 
the DMRG method, as described in Section IV E. The 
DMRG results do not give a principal advantage over 
QMC findings, but they provide a strong independent 
verification of the observed form of the spin gap, namely 
an exponential suppression at small J±. 



In Appendix A, we analyze the qualitative changes in 
the spectra of symmetric and single-pole ladders within 
a Jordan- Wigner mean-field calculation. Technicalities 
of spin-wave theory for long-range interaction are given 
in Appendix B. We finally present our conclusions in 
Section V. 
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II. PROBLEM SETUP AND QUALITATIVE 
CONSIDERATIONS 



We study the low-energy physics of the SSHL system, 
characterized by the following Hamiltonian: 

H = J||^(S M -Si, i+ i+cos 2 (§)S 2l i-S 2 , i+ i) 

i 

- J±^S M -S 2;i . (1) 

i 

Here, S Qji is a spin-1/2 operator acting on leg a and 
lattice site i and Jii > sets the energy scale. The single- 
pole ladder model 20 ' 23,25 corresponds to 9 = tt. 

It is convenient to reformulate the model Eq.(l) in 
terms of new variables, 

Si = Si,i + S 2 ,i , Ri = Si,i — S 2i j , (2) 

defined on the rung. The Hamiltonian then reads 

U = i^(l + cos 2 (f))(S J -S t+1 +R 4 .R t+1 ) 

i 



J± 



E ( s < - R * ) 



(3) 



The set of operators Si and Ri fully defines the o 4 
algebra. 20,21 

Let us define the retarded spin response function in the 
ladder situation 



Xji 



poo 



, Sji{q,uj') 
uj - ui' + iO 



(4) 



with j, Z = 1,2. At zero temperature the dynamic 
structure factor (spectral weight), given by Sji(q,u) = 
—ir~ 1 lmxji(q,u)) is represented as 

Sufaw) = J2(0\S?(q)\n)(n\S*(-q)\0)8(E n - E - u) 

n 

(5) 

where |0) stands for the ground state with the energy E , 
the sum runs over all eigenstates \n) of the Hamiltonian 
with energies E n and Sf (q) is a Fourier transform of Sf n . 
We also define the symmetrized combinations, 

S(q,uj) = Su(q,u) + S 22 {q,u) + Si 2 {q,uj) + S 2 i(q,ui), 
R{q,uj) = Su(q,u) + S , 22 (g,cj) - Si 2 {q,uj) - S 2 i(q,u>), 

(6) 

which are response functions for operators S z and R z in 
(2), respectively. 



Below we also use imaginary time correlation functions 



(7) 



(S z q (r)Sl q (0)) = J du i _ e _ f}u S{q,u) 
(R z (r)R z _ q (0)) = J du^^R&u,) 

with P = l/T. 

Let us now qualitatively discuss the situation. In the 
strong coupling region, J±/J\\ — > oo, and for all possi- 
ble twist angles 9 triplets on the rungs become more and 
more favorable such that the Hamiltonian of Eq. (3) re- 
duces to a pure spin-1 effective Hamiltonian: 



%eff — ^eff E] ^ i ' 

i 

Jeff - i(l + cos 2 ((9/2)) . 



(8) 



In units of effective coupling J e g (and thus irrespective 
of the twist angle) the spin gap of our model scales to the 
Haldane gap A H /J eS = 0.41048(6). 16 

In the weak coupling region the situation is more deli- 
cate and depends on the twist 9. Let us discuss here the 
limiting cases. 

For the symmetric ladder, 6 = 0, the gap opens as 
A <~ Ji, as obtained by bosonization and also quali- 
tatively reproduced in the simplified mean-field picture 
below. For the fully asymmetric single-pole ladder 9 = tt 
the mean-field calculation predicts a sub-linear behav- 
ior of the gap A ~ J±/J\\- The bosonization treatment 
becomes problematic in this case, as we cannot apply 
the continuum approach to the chain of dangling spins 
attached to the main leg. The assumption of the finite 
Fermi (sound) velocity in the subsystem breaks down and 
Jj_ cannot be used as a perturbation in the implicitly as- 
sumed hierarchy of scales J± <SC J\\ cos 2 (#/2) < Jy. 

Instead, we should start with the picture of a degener- 
ate band of dangling spins, whose degeneracy is lifted by 
indirect exchange between these spins through the main 
leg. This phenomenon is known cLS cL Suhl-Nakamura 
interaction 26-28 for the case of a dilute system of extra 
spins in a magnetic host. It has a direct analogy with 
the RKKY interaction where itinerant electrons medi- 
ate a long range spin-spin interactions between localized 
spins (see Ref . 29 and references therein) . In second order 
perturbation theory the SN interaction reads: 



J S n oc J 2 x(<7,w = 0) 



(9) 



where x(g, w = 0) is the spin susceptibility of the spin- 
1 /2 Heisenberg chain. We thus arrive at a similar energy 
J±/J\\ but now it stands for the SN- induced bandwidth, 
not for the spin gap. 

It is important to remark that both the mean-field 
treatment and the SN energy scale estimate provide us 
only with an upper bound for the gap value, J±/J\\- 
However this is likely a strong overestimation as quan- 
tum fluctuations, which are important in ID, are not 
accounted for. In Sec. Ill we present analytic arguments 
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in favor of a gap vanishing faster than with a power law. 
In Sec. IV we show, that the whole body of numerical 
data supports this scenario. 

III. ANALYTIC APPROACHES AND 
INTERPRETATION 

In most cases, it is instructive to map the system of 
spins onto a system of ID spinless fermions. We show 
in Appendix A, that such Jordan- Wigner transforma- 
tion and a subsequent mean-field theory analysis predict 
qualitative changes in the behavior of the gap with 9. 
Particularly, Eq. (A9) there indicates that the prefactor 
in the linear law, A ~ Jj_, diminishes with decreasing 
cos 2 , and that in the extreme case of single-pole ladder 
(cos | — > 0), the gap vanishes faster than | Jj_|. These ob- 
servations are qualitatively confirmed by our numerical 
simulations. At the same time, the attempt to compare 
the mean-field results for the fermionic theory with the 
gap extracted from QMC data shows that the mean-field 
theory overestimates the gap by an order of magnitude. 

In view of this fact, we perform a different type of 
analysis in the next subsections. 

A. Decimated blocks and effective spins 

Below we propose a scenario, which assumes different 
behavior of spin dynamics at high and low energies, as 
separated by the scale J± . This scenario is explicitly for- 
mulated for the single-pole ladder, which represents the 
most intriguing and difficult case, as seen in the numer- 
ics. The extension of this scenario to 9 ^ ir is also briefly 
discussed at the end of this subsection. The proposed 
theory allows to semi-quantitatively explain all features 
observed in the large-scale QMC and DMRG studies. 

For the high energies in consideration e > J± , the dan- 
gling spins should be considered as freely attached to the 
main chain. These energies correspond to short times, 
t < JJ 1 , and distances, x < J^/Jj_. Alternatively, one 
may think in terms of a higher temperature T > J±, 
which smears all fine features of the spectrum and leads 
to exponential decay of correlations beyond the temper- 
ature correlation length ~ J\\/T- The inverse length 
scale along the leg, associated with the crossover 

to the low-energy dynamics is defined by the relation 
J± ~ showing that soft spinon excitations in the 

HAF spin-1/2 chain, with momenta q < £ _1 are strongly 
intertwined with triplet-singlet transitions on the rungs. 

At these shorter time scales the classification in terms 
of singlet and triplet on the rung is not very appropriate. 
It means that the dangling spins are viewed as almost 
decoupled from the main chain: the situation is best de- 
scribed in terms of dangling spins coupled to each other 
by the Suhl-Nakamura interaction. The long-range char- 
acter of SN interaction leads to an almost flat dispersion 
of the magnon excitations in the subsystem of dangling 



spins, the estimated SN energy scale not exceeding the 
crossover energy Jj_. 

At smaller energies, e < J±, a picture of already 
formed triplets on the rungs is more appropriate. The 
discussion of the dynamics reduces to the rotations of 
the effective spins 5=1. Further, these rung triplets 
are interlaced by the leg interaction into large effective 
spin blocks of size £. Despite a possibly large spatial size 
of the block, the AF character of the main leg interac- 
tion selects the smallest possible total spin state of this 
block. Discarding the non-magnetic singlet, we focus on 
the block triplet state, i.e., when £ spins 5=1 are com- 
bined into a new effective spin 5 = 1 of the block. As a 
result, we have a model with nearest-neighbors interac- 
tion between large blocks. Such a model should display a 
Haldane gap, whose value can be determined from usual 
arguments. 

Let us start with the shorter distances and higher en- 
ergies. In this case the SN interaction between dangling 
spins 52, i and S2,i+ X has the form V(x)S2,iS2,i+ x with 

v(x) = Jlj ^x(q," = oy qx ~ Jl/M-iTHt/ x ) 

^ (10) 

with x(q,ui) is the response function for the HAF 5 = 
1/2 model. 6 For our purposes it suffices to approximate 
x(q, co = 0) by J^ 1 | cos(g/2)| _1 . The 1/q singularity near 
q = 7r shows that V(x) is sign- reversal and logarithmi- 
cally decaying with distance, V(x) ~ (— l) x ln(£/x). The 
scale £ in the argument of the logarithm, appears here 
as a parameter which will be further determined by a 
self-consistency criterion. Technically it is assumed that 
the allowed energies, co q , of the spin- wave continuum are 
restricted from below, cj q > Ju /£. 

A chain of dangling spins, coupled by long-range SN 
interactions (10), behaves differently from the standard 
HAF model with nearest-neighbor interaction. A simple 
spin-wave analysis is then indispensable here. Such an 
analysis cannot give a correct form of correlation func- 
tions, but delivers a qualitative information about the 
spectrum. 30 ' 31 

Using the formulas listed in Appendix B we conclude 
that the spin-wave spectrum in the system of dangling 
spins is given by the expression ujk — y/gk9k+v, with 
9k = V(ir) —V(n + k). Approximating the range function 
V(ir + k) ~ (j2/j || )( sin 2 (fc/2) +r 2 ) _1/2 we have 

Lo k+1T ~{Jl/J\\ )Zy] 1 - (1 + e sin 2 (fc/2))-V2 

|fc|»r 1 ( n ) 

Next we make the natural assumption that the top of 
the SN-induced band coincides (at least by the order of 
magnitude) with the logarithmic low energy cutoff intro- 
duced above, leading to {J±/J\\)£ ~ J\\/(, or 

S~J\\/J± (12) 
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Remarkably, this estimate shows that the low-energy 
spin- wave dispersion (11) is characterized by the same 
velocity LOk+ir ~ J\\\k\ as the higher-energy excitations 
in the main leg. This is despite the fact that, strictly 
speaking, the low energies e -C J± should be considered 
with a different approach as proposed below. Notice that 
the described spectrum resembles a simple picture of hy- 
bridization between the linear spectrum and initially de- 
generate band at non-zero energy <~ J± , with the crossing 
level repulsion phenomenon. 

At small energies e <C Jj_, we expect that the picture 
of spinons (in the main leg) scattering on the dangling 
spins, becomes inadequate. The dynamics of spins on 
the rung is characterized in terms of soft triplet dynam- 
ics, while the transitions to singlet state with the (now) 
large energy J± are discarded. On the same ground, 
we discard spin-wave continuum excitations with ener- 
gies e > J± <~ J/£, i.e., the description now has changed 
from the individual sites along the leg to entire blocks of 
length £. We may thus think in terms of the effective spin 
S = 1 on the rung, and these individual spins S = 1 are 
antiferromagnetically coupled to each other in the large 
block. 

Furthermore, we can characterize the whole block £ of 
spins 1 by its lowest non-trivial multiplct, a spin 1 again. 
However now it is an object defined at a much larger 
spatial scale. The trivial low-energy multiplct in such £- 
block is a spin-singlet state, which obviously drops out 
from the soft spin dynamics. 

The typical energy spacing between the resulting spin- 
triplet of the £-block and the higher multiplets is esti- 
mated again as Jy/£ ~ J±. We expect that this low- 
est triplet state is non-degenerate, i.e. other triplets are 
higher in energy. Below we refer to this lowest triplet 
state as ^-triplet. 

We note in passing that the estimate £ ~ J\\/J± fol- 
lows also from the argument presented in Ref. 22. It was 
suggested there that the spinons, propagating along the 
main leg by distance m and characterized by a typical 
energy J|| , break the pre-formed rung triplets resulting 
in an energy cost <~ mJ±. The resulting confining poten- 
tial should, in principle, restrict the motion of spinons to 
distances ~ J\\/J±- 

Our way of constructing ^-triplet resembles Kadanoff 's 
decimation procedure in the description of critical 
phenomena. 32 The new lattice of large blocks contains 
spins 1, which are denoted below by S^ n (here n num- 
bers a position in a new lattice) and are coupled to each 
other by a nearest- neighbor interaction. In fact, only the 
edge spins of each £ block are responsible for this inter- 
action. Adopting the notation that the weight of these 
edge spins in the ^-triplet is w\ <C 1 at £ ^> 1 (see below), 
we write explicitly for odd £ 



j+n 



(13) 



where n = [_j/£J and |_- ■ - J stands for the floor function. 
Similarly, = w 2 S£ n (-lp' +n . The phase (-l)J ac- 
counts for the AF character of the contributing spins in 



the £-block and the additional phase shift (—1)" is intro- 
duced to restore the translational invariance in the blocks 
picture. 

The definition of the above weight Wj is as follows. As- 
sume that the lowest multiplet in the £-block is a triplet 
\T,m), spanned by the spin-1 operators . Then, up 
to a sign, the edge operators of the block act as 

(T,m'\Sf tj \T,m)= Wl (T,m'\Sl n \T,m), (14) 

and similarly for the weight w 2 of S%j . The block Hamil- 
tonian on the decimated lattice reads 



.n'S'f ,n+l 



(15) 



with n = 1, . . . This model should exhibit a Haldanc 
gap, but the value of this gap is strongly diminished. 

Employing a standard albeit simplified approach used 
in the original Haldane paper, we apply the linearized 
spin- wave theory outlined in Appendix B and obtain the 
magnon dispersion, Wfe = Aw\ Jm sin(fc), where k is the 
wave vector on the decimated lattice k = 2im(^/L) with 
n = 1,2, It implies that the velocity of the lowest- 
lying excitations with respect to the real lattice is given 
by v — 4w 2 £J||. With the plausible assumption that this 
low-energy estimate coincides with the high-energy one, 
we obtain v ~ Jii and hence w\ ~ Notice that it 

corresponds to a bandwidth 4w 2 J|| <~ Jj_, in accordance 
with Eq. (11). 

We require that the zero-point magnon fluctuations 
cancel the local magnetization in the usual formula Sf j — 
1/2 - SijSfj, i.e., we write (S^-S^.) = 1/2. For the 
relevant low-lying modes, we have S^j = WiV^d^, where 
«4 is magnon creation operator in the nth ^-triplet. We 
then obtain the relation 



1/2 ~ 2w 



- ' 2 dk 

IT 



sinfc 



- 1 



where we introduced the cutoff wave- vector q , related 
to the Haldane gap A = 4u> 2 J|| sin<jo- This leads to the 
estimate qo <~ cxp(— 7r/4w 2 ) and 



A <~ w\ Jm exp 



(16) 



Let us discuss the decay of correlations at the distances 
r >• £ ; we take r/^ = x integer for simplicity. We have 

= {-l) r wl [ V — l- cosk „i(*+fc)» 
J -« 4?F ^ 2 + sin 2 fc 



7T 



(17) 
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The above assumption that w\ <~ £ 1 leads to the 
following final formulas 



(18) 



with C2 ~ 1. Our scenario also suggests that the long- 
range behavior of correlations (18) takes place for both 
main leg and dangling spins at distances r > £ ~ J\\/J±- 

Comparing the above formulas (18) to our numerical 
findings below, see Fig. 6 and Fig. 7, we verify that, in- 
deed, Wj oc J±. 

Finally, let us briefly comment on the situation where 
there is a weak exchange along the second leg, 8 ~ ir and 
J2 = J\\ cos 2 I -C J11 . Repeating the steps of the analysis 
leading to Eq. (11), we observe that the SN induced band 
appears on the background of usual exchange. Roughly, 
one can write V(k) ~ — J2 cos k + J]_/ Jy | cos(A;/2)| _1 and 
the bandwidth is estimated as J2 + £,J]^ / J\\ ■ This scale 
should not exceed the low-energy cutoff for the otherwise 
long-range SN interaction, which leads to the cor- 

rected estimate, £ _1 > max[Jj_/ J|| , cos 2 |]. The above 
scenario remains valid, as long as £ _1 <C 1, which im- 
poses restrictions on 8. We see in Fig. 8a, that the gap 
at 8 = 87r/9 (cos 2 | ~ 0.03) behaves qualitatively similar 
to the one at 8 — n, whereas the behavior of the gap 
at 9 = 7it/9 (cos 2 | ~ 0.12) is apparently different and 
closer to the one of the symmetric ladder, 8 = 0. For 
our semi-quantitative level of consideration these conclu- 
sions appear consistent and satisfactory. It was implied 
in Ref. 23 that the scaling of the gap with J± changes at 
a critical 8 C ~ tt. The above consideration suggests that 
there is no critical 8 but rather a smooth crossover be- 
tween two regimes, cf. also the estimate 8 ~ 0.547T in the 
next subsection for the crossover in strong rung coupling 
behavior. 



B. Strong rung coupling limit 

Let us consider the limit of strong rung coupling, 
|J_i_| 3> J\\- In the zeroth order of the small parame- 
ter J|i/| J_|_|, we have the effective Hamiltonian, (8). The 
perturbation is given by the operators Rf in (3), which 
connect the spin-1 rung sector, spanned by operators Sf , 
to the high-energy singlet rung state. The perturbing 
part is given by 



Hi: 



E 



(JrqR"Q" + JrrR°;R" + t) 



Q? — (Sf-i + ^i+i) 
M = i^||(l-cos 2 (f)), 



(19) 



and Jrr = + cos 2 (6>/2) = J e g. The leading correc- 
tion to the effective Hamiltonian (8) is obtained in second 



order of perturbation in Hint, by considering the virtual 
transitions to singlet states separated by the energy |Jj_|. 
This derivation is similar to the one of the t-J model from 
the large-?/ Hubbard model. We have: 

_n/(lA) ,nj(lB) 
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Using the identity 

R'jRj + SfSj = <W + itafi-fS'j 

one can eventually arrive to 
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The effective Hamiltonian (22) is rather complicated, 
containing quartic combinations of spins, and we propose 
here its qualitative analysis. The main term, Eq. (8), re- 
sults in antiferromagnetic correlations of the spins S{ at 
adjacent sites, whereas the next-to-nearest neighboring 
spins are aligned ferromagnetically. We thus expect Qi to 
be in the state Q = 2, and QiQi = Q(Q + 1) = 6. Next, 
let Qi and Sj be added into a multiplet of total spin p, 
so that S 4 Q, - § ((Sj + Q,) 2 - S 2 - Q 2 ) = \p{p + 1) - 4. 
One can check, that [QiQ, — S^Qi — (S^Q^) 2 ] = 0,6,0 
for states with p = 1,2,3, respectively. It means that 
the mostly AF correlated state with p = 1 obtains a 
higher energy due to i-L^\ while the less antiferromag- 
netically aligned state p — 2 leads to an energy gain 
— 6Jrq/|Jj_|. Alternatively, we may use the second line 

in the representation of ti^f^ in (22) and estimate it as ~ 
2 [x + 2x 2 - 1] J% Q /\J±\. Here S,S i+ i = x = -2, -1, 1 
for the total spin of the pair, Si + Sj+i, equalling 0, 1, 2, 
respectively; we also approximate Si_iSi+i ~ 1. Com- 
bining it with V. cS , we have the estimated energy per 
unit cell 



2 J 2 

SE ~ Tt\ 
\ J ±-\ 



Y(x + 2x 2 ) - 
W2)V 



(23) 



1 + cos 2 ((9/2) J 



In a simplified picture, we can compare the energy dif- 
ference, AE t8 , between the bond triplet state, x = — 1, 
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and the bond singlet state, x = —2. From (23) it follows 
that this energy difference per unit cell is 



(a) 



2?T f 1 27T 



AE, 



JeS 



2 J 2 
\J±\ 



-5F 



(24) 



For symmetric ladder, 9 = and Y = 0, the correction 
(23) favors the bond singlet and AE ts is always positive. 
For the single-pole ladder, 9 = it and Y = 1, the ex- 
pression (24) shows that the corrections ~ J 2 S /\J±\ are 
important even for Jj_ ~ 10J c ff, due to a large numerical 
prefactor. The sign of the correction in (24) changes at 
Y = 3/20 or 9 ~ 0.547T, and indeed we observe a slower 
saturation of the gap at J± 3> Jn for 9 > ir/2. 

Roughly, we can regard AE ts as a new value of J e g in 
Eq. (22), and it implies that the gap at 9 > tt/2 should 
follow the law A ~ 0AU eS (l-c(6)J e g/\ J±\) with c(9) ~ 
1. At large value of c(9), the intermediate region 1 < 
J±/J\\ S c (^) becomes rather extended. The detailed 
description of A(Jj_) at intermediate J± is beyond the 
scope of this study. 



IV. NUMERICAL ANALYSIS 
A. Exact Diagonalization 

In this section we analyze the SSHL model by means 
of exact diagonalization (ED) methods using the Lanc- 
zos algorithm. 33 ^ 34 Even though ED methods are limited 
to small systems, they provide considerable insight. We 
start our analysis with a study of the spin excitations 
Sji(q,u), Eq. (5) 

Fig. 2 presents the spin excitation spectrum for the 
isotropic ladder (9 = 0) and the single-pole ladder model 
(9 = 7r) in the weak coupling region. Precisely, it shows 
the dynamical spin structure factor (depending on the 
momentum q along the ladder), separately for the bond- 
ing (S(q,ui)) and anti-bonding (R(q, ui)) configuration. 
For the ladder, the dynamical spin structure factor for 
both bonding and anti-bonding cases displays features of 
the two spinons continuum of a single spin-1/2 chain 35 . 

Below we show, that such a continuum is qualitatively 
well reproduced by a mean-field theory, and corresponds 
to the particle-hole continuum stemming from the effec- 
tive fermionic Hamiltonian. The continuum for bonding 
combination, S(q,ui), is characterized by a slightly lower 
energy due to the weak ferromagnetic coupling between 
the chains. 

As apparent from Fig 2b, a narrow band emerges for 
the single-pole ladder model. We define its width, W, as 
the differences in energy between the low energy maxima 
at q = ir and q = tt/2 in the spin excitation spectrum 
for 9 = it. Associating this band with the SN splitting of 
dangling spins, we expect the width, W, to scale as J]_/J\\ 
in the weak coupling region. This is indeed confirmed for 
small system size in Fig. 3 where the ED data is found 
to fit well to a W oc J\ form. 



(b) 




FIG. 2: (Color online) Dynamical spin susceptibility in the 
weak coupling limit at J±/J\\ = 0.1: (a) ladder system (9 = 
0); (b) single-pole ladder system (8 = n). In both cases, 
results are obtained on 2 x 12 lattices with ED techniques. 
We choose a broadening s = O.Uii. The red (solid) lines 
represent the bonding spectrum, S(q,u), the blue (dashed) 
lines corresponds to the anti-bonding spectrum, R(q,uj). The 
spectral functions in the left panels are normalized to the 
structure factors S(q,t = 0),R(q,t = 0), respectively; the 
latter are shown in the right panels. 



For reasons that will be clarified below, we also inves- 
tigated the zero temperature spectral functions Sii(q,co) 
and 522(9, w), Eq. (5), with the emphasis on the weight, 
or residue, of the lowest excitation energy. This weight is 
a measure of the overlap between Sf (q)\0) and the first 
low lying excitation as modeled with the effective Hamil- 
tonian of Eq. (8). Fig. 4 plots this quantity, normalized 
by r o °° dwS u (q,u), that is: 



l(i|^(g)|Q)l 2 

(Sf(-q)S?(q))> 



(25) 



where the state |1) corresponds to the first magnetic 
excitation at wave number q. As apparent from the 
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FIG. 3: (Color online) ED results for the width W in the 
single-pole ladder model (0 = 7r). The effective SN interaction 
yields a bandwidth proportional to J±/J\\. 

ED results for the single-pole ladder model (see Fig. 4), 
Z^{q = it) corresponding to the second leg is almost in- 
dependent of Jj_. Similar results are found for other mo- 
menta. It means that nearly all the spectral weight for 
spins in the second leg belongs to the well -defined lowest- 
lying excitation. This should be contrasted with the sit- 
uation of spins the first leg, where only a fraction of the 
spectral weight belongs to this model, the rest partici- 
pating the spin-wave continuum. Clearly, the continuum 
fraction in Sn(q,u)) should increase as Jj_ is decreased 
up to the decoupled situation J± = 0, where one expects 
Z 2 (q = 7T)=0. 
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FIG. 4: (Color online) Normalized residue as a function of 
J±/J\\ for the ferromagnetic single-pole ladder. The calcula- 
tions were carried out with ED on an 2 x 8 lattice. 

We have also computed the spin gap A as a function 
of the interleg coupling J± for different values of 9. Un- 
fortunately, the finite size scaling becomes difficult in the 
weak coupling limit for all 9 (data not shown) and an 
extrapolation to the thermodynamic limit is not feasible. 
We can only confirm the differences in the scaling be- 
havior of the spin gap between the the single-pole model 
and the ladder model, where it is widely accepted that 
the gap opens linearly with the coupling between chains 

4,36 

Concluding this subsection, we notice that the advan- 
tage of the ED method is its high accuracy, but the 
method is limited to relatively small system sizes. For 



the single pole ladder, when the gap becomes small and 
comparable to interlevel spacing ~ J|| / L, the accuracy of 
the calculation becomes irrelevant. 



B. Quantum Monte Carlo, spin correlations 

To extend our analysis to larger system sizes, we also 
used quantum Monte Carlo (QMC) methods, performing 
simulations at finite inverse temperature (3 = 1/T. We 
applied two variants of the loop algorithm. For the spin- 
spin correlation function and for the string order parame- 
ter, discussed below, we used a discrete time algorithm. 37 
From the spin-spin correlation function we can then ex- 
tract the spectral function via stochastic analytical con- 
tinuation schemes. 38,39 In the next subsection, we also 
use a continuous time loop algorithm to directly com- 
pute the spin gap. 

We start our QMC analysis with a discussion of the 
dynamical spin-spin correlations, which can be computed 
on much larger sizes than with ED. However, the energy 
resolution is limited and hence we can only use this ap- 
proach at larger couplings than the ones reached with 
ED. 

The QMC results of the dynamical spin susceptibility 
for 9 = (ladder) and 9 = tt (single-pole ladder) with 
2 x 100 sites at J±/J\\ = 1.0 (/3Jy = 200) are shown in 
Fig. 5 (these results were partly shown in Fig. 4 of Ref. 23 ). 
At 9 — 0, inversion symmetry Si,j o S 2 ,i is present, such 
that the bonding and antibonding combinations do not 
mix. Since is even under inversion symmetry (with 
respect to the transverse direction), S(q,w) picks up the 
dynamics of the triplet excitations across the rungs. For 
ferromagnetic rung couplings Jj_ > 0, the low energy 
spin dynamics of the model is apparent in S(q, u>) which 
in the strong coupling limit maps onto the spin structure 
factor of the Haldane chain. In contrast, is odd under 
inversion symmetry and picks up the singlet excitations 
across the rungs. As apparent from Fig. 5b those excita- 
tions are located at a higher energy scale set by J± in the 
strong coupling limit. For the single-pole ladder model, 
9 = 7r, a mixing of the bonding and anti-bonding sectors 
occurs. As apparent in Fig. 5c, d both R(q, lo) and S(q, ui) 
show high and low energy features. The low energy, nar- 
row, dispersion curve in Fig. 5c is a consequence of the 
SN interaction and reflects the slow dynamics of triplets 
formed across the rungs. 

In spite of the limited energy resolution of the QMC 
method, we can extract the value of the spin gap, by 
studying the spin correlation functions at long distances. 
Such analysis also provides an insight of the intermediate 
energy scales. We show the behavior of the correlations 
( S l,i S t,j) and ( S 2,i S 2,j) as function of \i - j\ in Fig. 6 
for a particular value Jj_ = 0.5Jy ; the distance is mea- 
sured in lattice spacings. We notice that beyond a certain 
length scale both correlation functions behave similarly, 
differing only in overall factor. We hence plot in Fig. 
6 the function (Sf^Sfj) as is, while (S^iS^j) is mul- 
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(a) S(g,u) for 9 = 




FIG. 5: (Color online) QMC results of the bonding and anti- 
bonding dynamical spin susceptibility for the ladder (6* = 0) 
(two top panels) and single-pole ladder (6 — it) (two bottom 
panels) systems at J±/J\\ = 1.0. /3Jii = 200 was taken in the 
simulations. 



tiplied by a factor discussed below; after this "normal- 
ization" both curves are indistinguishable at large dis- 
tances d > 10 (these results were partly shown in Fig. 3 
of Ref. 23 ). This indicates that the lowest-energy dynam- 
ics of dangling spins S| i and of the main leg S* i is the 
same slow dynamics of rung triplets, wherein these spins 
enter with different weight. 

At largest distances an exponential decrease of correla- 
tions is observed, corresponding to a gap in the spectral 
weight of Sjj(q,io). In the previous section we devel- 
oped a theory which accounts in a semi-quantitative way 
the whole body of numerical findings. According to this 
theory, the long distance behavior of the spin-spin corre- 
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FIG. 6: (Color online) Equal time spin-spin correlation func- 
tion at J±/J\\ = 0.5 for the single-pole ladder along both 
chains. Simulations were done at /3Jy = 5000 on a 2 x 800 
lattice. 

lations is given by 

w? /Ar\ 
\(SlrSlo)\ = i- K °{-)- ( 26 ) 

with Kq the modified Bessel function, v ~ Jn velocity of 
spin excitations, Wj (with j = 1, 2) is the weight of the 
spin Sj r in the effective spin-1 variable of the single-pole 
ladder. As explained above, we expect wj ~ J±/J\\ while 
A ~ J L exp(-Jy/Ji). 

In Fig. 6 we fit the long-ranged equal time spin-spin 
correlation function at Jj_/J\\ = 0.5 to this form. Several 
comments are in order. 

a) Normalizing the spin-spin correlations of the dan- 
gling spin by a factor 2.18 = w\jw\ provides a per- 
fect agreement between the long range correlations on 
both legs. Note that the numerical factor 2.18 is close to 
^2(3 = n)/Zi{q — 7r) ~ 2.25 as obtained from the data 
of Fig. 4. Hence the low energy dynamics of the spins on 
both legs are locked in together. This observation con- 
firms the picture that the low lying spin mode observed 
in Fig. 5c indeed corresponds to the dynamics of triplets 
across the rungs. 

b) One can read off a length scale at which the func- 
tional form of both correlation functions differs. This 
length scale marks the crossover from high to low energy 
beyond which an effective low energy theory can be ap- 
plied. 

c) In our previous paper 23 we also pointed out the exis- 
tence of the intermediate asymptotic, \i — j| -1 ^ 3 , in case 
Jx = 0.3 Jn . The existence of such power-law decay is in- 
teresting on its own, but the detailed description of this 
regime is beyond the scope of this study. This is the case 
when the gap becomes comparable to the energy spacing 
due to the finite size of the system, and it causes diffi- 
culties in other complementary numerical techniques, see 
e.g. DMRG results below. 

d) For comparison we have plotted the spin-spin correla- 
tions of the spin-1/2 chain 22,40 which, on a length scale 
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set by the correlation length, decay much faster than the 
correlations of the single-pole ladder. This very slow de- 
cay should be seen as a direct consequence of the long 
ranged nature of the SN interaction, 
c) The fit to the form of Eq. (26) is next to perfect 
thereby providing an excellent description of the low en- 
ergy physics. The ratio of the gap to the velocity as 
well as the weight of the spins in effective spin-1 vari- 
able are plotted in Fig. 7. In particular, assuming a 
linear dependence of the weight (see below), we obtain 
w\ = 0.13J_i_/J|i and w\ — 0.060Jj_/J|i- The gap values 
are fitted by an exponential law, A ~ Jx exp(— Ju/Jx), 
as discussed in Sec. Ill A. In the next subsection, we will 
see that this exponential form provides an excellent fit to 
the spin gap directly measured in QMC. The parameters 
of the fit shown by a dashed line in Fig. 7 are explained 
in the last paragraph of the next subsection. 



0.03 




0.2 0.4 

J±/J\\ 
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FIG. 7: (Color online) Scales in the correlation function (26) 
on the second leg as a function of J± / Ju . The dotted line is a 
fit to a linear law and the dashed curve is a fit with v — 0.28 Jy . 
The form of A is discussed at the end of Sec. IV C. 



C. Quantum Monte Carlo, spin gap 

For the spin gap calculation the continuous-time loop 
algorithm of ALPS 41 was used. Here, we can calculate 
the correlation length in imaginary time £ r (?) for a given 
wave vector q via a second moment estimator 16,42 : 



Mo) 



2tt 



X(q,w = 0) 
X (<7, W = 27r//3) 



1/2 



(27) 



with x(o.i w ) = Jo dTe lLJT x(<li t) the Fourier-transform of 
the imaginary time dynamical structure factor (7). 

The inverse of £ T (q) converges in the limit L, (3 — > oo 
to A(q) = - E , if A(q) ^ 0. Here E x {q) is the 

minimum of the dispersion at wave vector q and Eg the 
ground-state energy: the spin gap is simply obtained as 
A = min g A(g). If the system is gapless at q, £ T (<7) _1 
is an upper bound of the finite-size gap at q for any fi- 
nite L and ft. The full dispersion curve can therefore 



be calculated in principle with this method. In prac- 
tice however, the simulations suffer from large statisti- 
cal errors in %(g, w) when q is different from the wave 
vector of the lowest lying excitation. The situation can 
be ameliorated by using improved estimators 43 for the 
imaginary time dynamical structure factor, which is sim- 
ply related to the loop sizes in the algorithm. The wave 
vector picked by the loops in the algorithm corresponds 
to the one given by the sign of the coupling constants 37 , 
which is in our case q — (tt, 0) (for ferromagnetic Jj_ > 
and antiferromagnetic Ju < 0). The improved estimators 
have a smaller variance than conventional ones, and we 
therefore obtain good statistics for £ r (7r,0) and the spin 
gap A = A(ir, 0) when it is finite. 

Now we turn our attention to the spin gap as a function 
of the coupling J± and twist angle 9 (see Fig. 8a) (these 
results were partly shown in Fig. 2 of Ref. 23 ). 

In the strong coupling limit Jx/Ju — > oo the model 
maps onto the AF spin-1 Heisenberg chain, Eq. (8) and 
we expect the Haldane gap Ah — 0.41J c ff . This behavior 
is clearly apparent in Fig. 8a. However, as 9 grows from 
9 = to 9 = tt the approach to the Haldane limit becomes 
slower. Whereas the scaling for 9 = and 9 = tt/2 
are nearly similar, the scaling behavior for the single- 
pole ladder (9 = tt) differs, as can be seen in Fig. 8c. 
In Section IIIB we argued that at largest Jx the gap 
should follow the law A ~ 0.41J off (l - c(0) J eS /\ J±\). 
Now we confirm this picture by fitting the numerical data 
with this form at the largest Jx, as shown in Fig. 8c. 
We see that c{9) definitely grows as a function of 9 and 
reaches a large value c(9) ~ 8 for the extreme case of 
the single-pole ladder, 9 = tt. In the intermediate region, 
1 < Jx/J\\ < c(0), we cannot expect a good fit of the 
gap, as clearly visible in Fig. 8c. 

For the ladder system (9 = 0) our data stand in agree- 
ment with the independent QMC calculations of Ref. 36; 
the spin gap opens linearly with respect to the coupling 
Jx up to logarithmic corrections. It is beyond the scope 
of this work to pin down the exact form of the logarithmic 
corrections, and we refer the reader to Ref. 36 for further 
discussions. This behavior is stable up to to large twist 
angles, and it is only in the very close vicinity of 9 = tt 
that a different behavior is observed. 

As discussed in the previous section, the spin gap at 
9 = tt is expected to decrease exponentially with de- 
creasing values of Jx, as oc J± exp(— Jii/Jl). In or- 
der to fit the QMC data of Fig. 8b in a wider region, 
Jx % J\\, we also allow for a correction in the pref- 
actor in this law, namely we assume the dependence 
A = aJx(l — bJx/J\\) exp(— cJu/Jx)- Fitting the data 
of Fig. 8b to this form gives a = 0.077, b = 0.32, c = 1.34. 
These parameters can now be used to fit the data for 
the gap, A/v, as extracted from the spatial decay of 
correlations in Fig. 7. The only adjustable parameter 
there is v/Jn, and we obtain a good agreement in Fig. 
7 for u = 0.28 Jy. The effective model of the previ- 
ous section leads to an expression v = 4iu 2 £J||. Com- 
paring these values, we determine the crossover scale 
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FIG. 8: (Color online) (a) Size and temperature converged 
values of the spin gap A as a function of the coupling J± / Jm 
for various twist angles. The gap is rescaled by J c ff = 
4* (X + cos 2 (#/2)) such that in the large- Jx-limit it converges 
asymptotically toward the Haldane gap of a AF spin-1 chain 
A_sr / Jeff ~ 0.41. (b) Spin gap for the single-pole ladder model, 
8 — tt. Lines denote quadratic and exponential fits for the 
spin gap (see text), (c) Spin-gap for SSHL at larger J± for 
different 8s, see detailed discussion in Sec. IIIB. 



£ = 1.2J|i/Jj_, which separates the long-distance behav- 
ior from the short-distance one. For the particular value 
J± = 0.5 Jy , we obtain £ ~ 3, and this scale is clearly visi- 
ble in the merging of curves for the spatial dependence of 
correlations in the first and second legs, as shown in Fig. 
6. We further show in Fig. 8b the result of a quadratic fit 
to the gap A ~ (J±/J^) 2 . While such a fit is reasonable 
for low J_l as remarked in Ref. 23, the above exponen- 
tial form appears to provide a better fit in a larger J± 
window. 



D. String Order Parameter 

To pin down the nature of the ground state, and in par- 
ticular for the single-pole ladder, we compute the string 
order parameter characterizing the Haldane phase. In the 
strong coupling region the system maps onto an effective 
spin-1 chain, for all twist angles 9 and the ground state 
can viewed in terms a valence bond solid (VBS) 5:44,45 . 
In the VBS state the spins on a rung form triplets in 
such a way that neglecting triplets with z-component of 
spin m = reveals a Neel order. This hidden AF order, 
is characteristic of the Haldane phase and, as shown by 
Nijs and Rommelse 46 , is revealed by the non-local string 
order parameter 



a = (^ exp 



no+L/2- 

• E 

j=n 



Si 



s. 



no+L/2 



(28) 



where Sf = S I i + 51 i , rig stands for an arbitrary rung 
and L denotes the system length. O s is also sensitive to 
a true AF order. To distinguish between a hidden AF 
order and a true Neel order another order parameter has 
to be introduced 46 



exp 



no+L/2 

■ E ss 

j=n 



(29) 



which is zero in the Haldane phase (hidden AF order) and 
finite in the Neel phase. Starting from the strong cou- 
pling region where the system is clearly in the Haldane 
phase, O s ^ and Oh = 0, we analyze the evolution of 
the string order parameter as a function of the coupling 
J± and the twist angle 6. For 6 = the order parameter 
O s stays finite and Oh is zero for all couplings. Hence 
the ladder system remains in the Haldane phase, inde- 
pendent of Ji. 

The situation for the single-pole ladder is much more 
delicate (see Fig. 9a, these results were shown in Fig. 5 of 
Ref. 23 ). For J±/J\\ < 0.4 the string order parameter Oh 
appears to be non-zero, thus indicating Neel order. In the 
previous section we have shown, that at weak couplings 
and for 9 = tt the SN interaction generates a very slow 
decay of the spin correlations. For instance, the correla- 
tion length £ for J±/J\\ = 0.3 is larger than the system 
size and the very slow decay of the spin-spin correlation 
functions at distances smaller than £ mimics AF order. 
However, when increasing system size, one expects Oh to 
vanish and O s to converge to a finite value. This expec- 
tation is supported by a finite size scaling of the string 
order parameters, as shown in Fig. 9b for 9 — 8tt/9 and 
6 = tt. The crossover between the AF order for small 
systems and the disordered phase is rather obvious for 
9 = 8tt/9 (cos 2 | ~ 0.03) at J x /J\\ = 0.2. For small 
lattice sizes the system seems to be in the AF ordered 
phase as indicated by the fact that both string order pa- 
rameters are finite. However, for increasing lattice sizes 
the order O s remains nearly constant, whereas the order 
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FIG. 9: (Color online) (a) String order parameters O a and 
Oh as a function of interleg coupling J± for 9 — n. In the pa- 
rameter range J±/J\\ < 0.5 finite size effects are still present, 
(b) Finite size scaling of the order parameters for 9 = 8n/9 
(J±/Ju = 0.2) and 9 = n (Jj_/Jn = 0.3). The simulations are 
carried out up to /3Jy = 7000 and 2 x 800 spins. 



parameter Oh decreases and finally vanishes. At this 
point the order parameter O s increases again. For the 
single-pole ladder (9 = ir) at J±/J\\ = 0.3 we do not ob- 
serve the disappearance of Oh, which shows that finite 
size effects are strong in this case even for a system as 
large as L = 800. 



E. DMRG analysis 

In this section we present our analysis of the Hamilto- 
nian, Eq. (1), with the use of the Density Matrix Renor- 
malization Group (DMRG) 47 for the single-pole ladder, 
i.e. — 7r. Overall, our results presented in Fig. 10 clearly 



support a non-analytic exponential scaling of the gap in 
Jj_ as suggested by the analytical considerations in this 
paper. 

For the calculation of the spin gap the specific choice 
of the boundary condition (BC) is crucial. While DMRG 
prefers open BCs for numerical stability and accuracy, 
they must be dealt with carefully in order to separate 
boundary effects from bulk effects. Periodic boundary 
conditions can be applied within DMRG, 48 - 49 yet with 
somewhat limited accuracy and efficiency. In this pa- 
per we therefore adhere to the conventional DMRG with 
open BC. In order to still deal with a Hamiltonian with 
periodic BC, a long bond connecting the ends of the chain 
can be introduced for short system sizes. Alternatively 
for somewhat longer system sizes, the chain with peri- 
odic boundary can be reshaped into a double chain with 
ends connected to form a loop. We adopted the latter 
approach since it is stable in finding the ground state of 
the system. Nevertheless it is enormously costly numeri- 
cally, and keeping up to 5120 states for chain lengths up 
to L = 256 rungs total, the DMRG results still showed 
significant uncertainties in the ground state energy for 
small J±/J\\ insufficient to accurately resolve the expo- 
nentially small gap (see Fig. 10, panel a). 

When compared to the other calculations (see, e.g., 
panel (d) below), the gap for periodic BC is consistently 
overestimated for small Jj_/J\\. This comes from the 
fact that the ground state is typically well-represented 
(smaller block entropy due to gap), while the excited 
state at larger 5* ot has a larger block entropy as it is 
a part of the continuum. Due to the limited number 
of states kept, the excited state is less accurately repre- 
sented which leads to an overestimated gap. This effect is 
more pronounced for large system sizes as can be clearly 
also seen from Fig. 10a. 

At the same time the numerical data for the gap for 
large J±/J\\ are reliable, allowing for an extrapolation 
towards the known value of the Haldane gap with less 
than 1% relative error as shown in panel (a). The fitting- 
range used was J\\/J± < 0.55 as indicated by the vertical 
dotted guide in the panel. With a fit of the form A( J±) ~ 
e -const J\\/Jx this shows that the Haldane gap is reached 
rather slowly when increasing J±. 

With periodic BC being of limited accuracy as ex- 
plained above, we adopted the plain open BC also on 
the level of the Hamiltonian for the rest of the DMRG 
calculations. Hence we deal with a single ladder, in con- 
trast to the connected double ladder above. Keeping up 
to 2560 states leads to clearly converged numerical data 
for all J±/J\\ analyzed. Note nevertheless, that the block 
entropy rises rapidly for Jj_/J\\ < 0.5 due to the near de- 
generacy of the dangling spins on the single-pole ladder 
in the limit J±/Ju — > 0. 



From the numerical analysis for systems with open BC, 



one observes spin 1/2 edge excitations visible in terms of 
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FIG. 10: (Color online) DMRG analysis of the single-pole ladder (8 = n) - Panel (a) shows the spin gap between the ground 
states of the S z = and S z = 1 spin sectors for several system sizes using periodic BCs (double chain, see text). The error 
bars indicate the convergence with respect to the number of states D kept (with an extrapolation in 1/D — > with D < 5120). 
The convergence to the Haldane gap for J± ^> Ju is reproduced within 1% relative error, within the fitting-range at small 
J\\/J± indicated by the vertical dotted line. Panel (b) shows spin gap data using open BCs for the ground states of several spin 
sectors S z £ {0, 1,2,3,4}. Keeping up to 2560 states, all energies are converged with negligible uncertainties. The consecutive 
energy differences (spin gap) A§ s^+i between the ground states of the spin sectors S z and S z + 1 are plotted vs. J±/J« for 
L = 128 (dashed) and L — 512 (solid), respectively. The vertical order (top to bottom) of the entries in the legend matches 
the order of the lines appearing in the panel (top to bottom). The extrapolated spin gap Ajf (grey and black dots for L — 128 
and L — 512, respectively), as well as A^ (red crosses), as obtained from a finite size scaling of Af j2 for systems of length 
L 6 {128,256,512}, are shown. Panel (c) shows Aq from panel (b) for L £ {128,512} and also for L = 256 together with 
A^j (all for open BC) on a semilog plot again with inverted x-axis similar to panel (a). The extrapolated value of the gap for 
small J\\/J± for open BC compares well with the exact Haldane gap within 1% relative error already for the smaller system 
size L — 128 (fit range up to first vertical dotted guide). Fits to the data for large J\\/J± are shown in red solid and dashed 
curves. The fit range chosen is again indicated by the two vertical dotted guides, here for Ju/Jx > 1. Panel (d) summarizes 
the data obtained with periodic (green solid) and open (red asterisk) boundary conditions DMRG calculations, as compared 
with QMC data (black solid). The red solid line is a fit for large J«/J± similar to the one in panel (c). Error-bars indicate the 
respective uncertainty in energy. Inset shows that A/J± reaches the maximum at J± ~ 2J\\. 



an alternating finite (S z ^) along the sites i which de- 
cay exponentially with the distance from the ends of 
the chain. 50 As the total spin Sl ot is a conserved quan- 
tum number, the individual symmetry sectors can be 
analyzed separately, with the overall ground state lying 
within 5' ot = 0. However, as it turns out, the Sl ot = 



sector has clear S z = ±| edge excitations with oppo- 
site signs being consistent with Sl ot = 0. The situa- 
tion is similar in the Sl ot = 1 sector, but this time yet 
with equal signs of alternation, consistently adding up 
to 5*°* = 1. Therefore 5*°* = and S"* ot = 1 are de- 
generate in the thermodynamic limit yielding a four-fold 
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degenerate ground-state due to the presence of the open 
boundary 50 . Only starting with 5*°* = 2 a true bulk ex- 
citation is generated. 51 . In order to extract the spin gap, 
we therefore prefer to monitor the splitting in energy with 
respect to S** * > 2. To this end the energy differences 
A| s +1 between the ground states of the consecutive 
spin sectors S z and S z + 1 are calculated and plotted in 
Fig. 10b. From the above argument, Aq 1 clearly vanishes 
in the thermodynamic limit as the overlap between the 
spin 1/2 excitations confined to the boundaries decays 
exponentially with system size. Therefore only A| s +1 
for S z > resembles the gap with a finite-size effect 
added with each increment of S z . Note in Fig. 10b that 
all curves for S z > start falling onto the same line as 
L increases. This indicates, consistently with the notion 
of a spin gap, that for large enough system sizes each 
increment of 5*°* by 1 costs an energy equal to the gap 
value. Note, however, that the lowest energy states for 
the S z > 1 sectors may already lie within a continuum of 
states. 

The strategy then to extract the spin gap is two-fold: 
(i) extrapolation of Af s +1 for S z > for constant 
system size L, using a quadratic fit towards S z = to 
eliminate finite-size effects with increasing S z , 



In ee lim Ao 



Sz+1 , (L = const), 



(30) 



(grey and black dots in Fig. 10b), and (ii) actual finite- 
size scaling on Af 2 , i.e. the lowest S z that yields a finite 
spin gap Af s +1 in the thermodynamic limit, 



lim Af 9 

1/L-vO ' 



(31) 



(red crosses Fig. 10b). As can be seen in the same panel, 
both strategies nicely agree with each other for L = 512, 
indicating that the data are well converged and consis- 
tent. It also shows that Eq. (30) is a valid way of ex- 
tracting the spin gap for systems that are large enough. 

In order to analyze the data at small J±/J\\, the data 
for Aq and A^ are (re)plotted in Fig. 10c with inverted 
x-axis on a semilog-y plot, as we expect a non-analytic 
behavior of the form A <~ e - const J i\/ J ± (and as motivated 
by previous studies on spin ladders too 52 ). As a consis- 
tency check for open BCs, the gap Aq for L = 128 is 
again extrapolated towards J\\/J± — > to retrieve the 
Haldane gap with reasonable accuracy of 1%. 

For large J\\/J± > 1, the gap A^, in the thermody- 
namic limit (red crosses in Fig. 10c) can be fitted nicely 
using exponential forms of either type 



A(J') = a e- Cl/J \ and 
A(J') = (oiJ' - a 2 (J'f) 



-C 2 /J' 



(32a) 
(32b) 



with J' = J±/J\\ , and ai and Cj being fitting parameters 
(see Fig. 10c where x = l/J'). Note that the (J') 2 term 
in Eq. (32b) is important to obtain a clear agreement 
with the numerical data.The data are equally well fitted 



by both forms in the region J\\/J±_ > 1, with deviations 
of the fit in Eq. (32b) (dashed line in Fig. 10c) visible 
only outside the fitted region, at Jn < J±. Note also 
that despite the fact that the fitting-range chosen to be 
J\\/J± € [1.05, 2.25] (as indicated by the vertical dotted 
guides in Fig. 10c), both fits extrapolate well up the last 
data point at 2.5. 

Finally, the DMRG results for the gap are summa- 
rized and directly compared to the QMC simulations in 
Fig. lOd, as a function of J±/J\\. The gap A^, in the 
thermodynamic limit (red asterisks) clearly lies within 
the error bars of the other less accurate data sets, while 
its own error bars are negligible. The results obtained 
this way are then reliable down to smaller J±/J\\. The 
exponential fit reproduced from Fig. 10c (solid red line) 
in panel (d), finally, illustrates the extremely fast decay 
of the gap towards small J± . Yet as clearly supported by 
Fig. 10c, the spin gap remains finite in this region. The 
inset in Fig. lOd illustrates that, when scaled to ,7^, the 
gap in the single-pole ladder has a maximum at J± ~ 2 Ju 
and decreases exponentially at smaller J±. This is to be 
contrasted with the symmetrical ladder where A/Jj_ sat- 
urates to a constant at J± — > 0. 



V. CONCLUSIONS 

In conclusion, we investigated asymmetric spin lad- 
ders, with different values of exchange interaction of spins 
5 = 1/2 along the two legs as parametrized by 9. For 
ferromagnetic rung coupling J± the spectrum of excita- 
tions is characterized by a Haldane gap, as expected for 
the effective spin-1 rung variables which are coupled an- 
tifcrromagnetically along the chains. We confirm this by 
the numerical analysis of the spin gap, spin correlation 
functions and of the corresponding string order parame- 
ter. 

The most intriguing behavior is observed near the 
single-pole situation, i.e. in the absence of exchange 
along the second leg. Our extensive numerical analysis 
shows that the spin gap decreases with Jj_ exponentially 
fast, A <~ J± cxp(— Jii/Jj_), unlike the conventional sym- 
metric ladder behavior, A ~ Jj_. In order to explain 
the whole body of numerical data, we develop a theory, 
which takes into account the indirect Suhl-Nakamura in- 
teraction between spins and the formation of large ef- 
fective blocks of spins. In a certain sense, the formula 
for the gap as obtained from this approach combines the 
"quantum" prefactor J± and the semi-classical exponent 
exp(— J\\/J±) arising from the large-block picture. 

In summary, we have a strong evidence that a spin gap 
opens for any J± but it may be exponentially suppressed 
nearly single-pole situation. This smallness leads to dif- 
ficulties in the actual observation of this gap, both due 
to the finite resolution of the tool employed (numerical 
or experimental) and to the finite-size effects associated 
with a large correlation length. 

The case of a negligible gap was reported in Rcf. 23 



15 



(a) (b) 

Leg 1 7r/2 tt/2 




FIG. 11: (Color online) (a) Zigzag path used for the Jordan- 
Wigner transformation, Eq. (Al), (b) gauge transformation, 
leading to Eq. (A3). 

for Jl = 0.3Jy. We did not discuss this situation in the 
present paper, since at these parameters the correlation 
length is larger than the system size and the true string 
order parameter is not formed, cf. Fig. 9. In such case 
the observed power-law decrease of correlation functions 
should be viewed as an intermediate asymptote, and the 
origin of the particular value of the decay exponent re- 
ported in 23 is not clear. 
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Appendix A: Jordan- Wigner Mean-Field Approach 

The definition of the Jordan- Wigner transformation for 
the ladder topology relies on the choice of a path labeling 
different sites. With the choice shown in Fig. 11a it reads: 



ceeding to the details of the mean-field calculations be- 
low, let us summarize the results obtained within this 
approach. The ground state of the asymmetric ladder 
is characterized by spinless fermions circulating around 
a plaquette thereby allowing for a 7r-flux phase as solu- 
tion of the mean-field equations. This leads to a spin gap 
A cx |JjJ cos(0/2) for |Jj_| < J\\ cos(0/2). Such regime is 
unavailable for the single-pole system at = it in which 
case we obtain an (indirect) spin gap A cx J\/J\\- The 
mean-field also predicts a smooth crossover between these 
two regimes. 

Our analysis is rather standard and similar to Ref. 53. 
After application of the Jordan- Wigner transformation, 
the Heisenberg Hamiltonian of Eq. (1) can be written as: 

8=1 

-4cos 2 (|)^((if) 2 -ife^), 

i=l 

i=l 

Here, we have defined 

A (Q) = c l,i c «,i+i + h.c. , B t = c f M c 2!l + h.c. (A2) 

We restrict ourselves below to a phase with zero magne- 
tization, (S^i) — 0, which corresponds to (n a ,i) = To 
proceed further we will replace n a ^ by its mean value. Al- 
though this simplification cannot be rigorously justified 54 
and more elaborate treatments are possible in a Jordan- 
Wigner approach 55 , we show below that it qualitatively 
reproduces the results available by more sophisticated 
methods. 

The replacement e mna - i+1 — > i defines a spiral Z7(l) 
phase shift for fermions on each leg. The relative phase 
between these spirals demarcate two qualitatively differ- 
ent situations. In one case the phase flux <!> through each 
plaquette is zero and in another case it is equal to it. Us- 
ing the gauge transformation, we can reduce the <I>-flux 
phase effective Hamiltonian to the form 



i-1 

S ti = c t lj4 exp(-i7r^(nij +n 2 ,/)) 
i=i 

i i—1 

St t i = 4,iexp(-i7r(^n 1)i +^n 2 ,i)) (Al) 
i=i i=i 

where n a ^ = ;C a j is the density operator at site i with 

leg index a = 1,2. i and c a ^ are spinless fermionic 
creation and annihilation operators which satisfy the an- 

ticommutation rules ■! c a> i, cjj ■ > = 5ij5 a p. Before pro- 



= E((a (1) ) 2 +a (1) ) 
i 

- J i^(l)j:((APr + AP), (A3) 
I 

i 

with a checkerboard character of coupling in the J± chan- 
nel happening in the 7r-flux case, $ = it. On the other 
hand, the 0-flux case is characterized by a uniform J± 
coupling, e 4 *' -> 1 in Eq. (A3) . 
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We then use the mean-field decoupling, 
{A ( f ) )=A^\ (B l )=e m B 
for the quadratic part of the Hamiltonian (A3), where 



d0')\3 



A\ i] sliould 



the special property Bf — Bi, (A 
be taken into account. The remaining Hamiltonian is 
quadratic in fermions and its spectrum is easily found. 
It can be shown that the 7r-flux phase provides a lower 
ground state energy, and hence we focus on this phase 
from now on. From the consistency equations below it 
can be shown that A^ = A^ = A, and we can use this 
observation to simplify our subsequent formulas. 
The spectrum has two bands, 



4 ±) = ^||(l + ^)sin 2 (|)cos g ± \E q 



(A4) 



E q = ^(J||(l + A)(l + cos 2 (f)) cos q) 2 + Jl(l + B f, 

and this dispersion should be combined with the consis- 
tency conditions 



A = 1 ^ J||(l + ^)(l + COS 2 (|))cOS 2 g ^ 



E„ 



B 



1 ^ J±(l + B) 



E 



2E n 



(A6) 



where the thermodynamic limit \ ^2 q — > / dq/2ir is as- 
sumed. 

In the limit Jx — > we have A ~ 2/ir and 



B ~ 



(Jx/J||)ln(J||/Jxci) 
(l + 7r/2)(l + cos 2 (f))' 



with ci <~ 1. 

At half-filling, corresponding to a vanishing total mag- 
netization in the spin language, a direct gap is given by 



e (+) - - 



(-) _ 



= £? 9 at g = 7r/2 or 

A =|J ± (1 + B)|, 

^|Jx|(l + 0((Jx/J||)ln(J||/Jx))) 



(A7) 



Hence, this simple mean-field approach is consistent (in- 
cluding logarithmic corrections) with bosonization 4 and 
quantum Monte Carlo simulations 36 for which a spin gap 
of the form Eq. (A7) is found in the weak coupling limit. 

The indirect gap, which is the minimum excitation en- 
ergy in the spin system at zero temperature, is defined 
as A = mm q s q + ^ — max,t e[ ' = 2xmn q e q + \ The qual- 
itative form of the spectrum (A4) at 9 ~ 7r is depicted 
in Fig. 2 of Ref. 21. For small J±, a fiat band is appar- 
ent reflecting the macroscopic degeneracy of the model 
at 9 = ir and J± = 0. This leads to a dense spectrum of 
particle-hole excitations at low energies. A straightfor- 
ward calculation shows that 



A a 2 cos I 
A = A 0: 



1 + cos 2 % 



(A8) 



at cos | ~ 1 and small enough J±. When cos | — > 0, 
the domain of linear dependence of A on Jj_ disappears 
and we have the quadratic law. Expressing energies in 
units of e = J||(l + + cos 2 |), we obtain for small 
Jx/eo, cos § <C 1 : 

A ~ 2 Jx cos f , Jx < 2e cos § 

~ Jx/Oo) + 2e cos 2 §, Jx > 2e cos f . (A9) 

A linear regime for the indirect gap occurs for incom- 
mensurate q in the above expression min g e q + ^ , whereas 
the quadratic regime in (A9) corresponds to the differ- 
ence e q t} n — max fc e k Z] , i.e. commensurate wave- vector 7r 
of the particle-hole excitation. 



Appendix B: magnons for long-range interaction 

We use the Dyson-Maleyev representation of spin op- 
erators. In a predominantly AF situation, we write 



Sf =(a - a}a«)(-l)' 

Sf =y/s/2(a\ + a; - ajaf/(2s)) 

S!=i^/2(al-a l + a\a?/(2s))(-iy 

For spins in different sublattices we have 

SiS 2 = -s 2 -s + s(a\ +a 2 )(4 + ai) + ... 

and for spins in the same sublattice : 



(Bl) 



SiS 3 



s(a\ - 4)(oi - a 3 ) + . . . 



Adopting the FM interaction V"(r) within one sublattice 
and AFM interaction V(r) between sublattices, we come 
to the linearized Hamiltonian : 

H = S ^4^[V(0)+V(0)-V(k)}+V(k)[alal k +h.c.}) 

k 

(B2) 

Notice that we have V(k + it) = —V(k) and V(k + ir) = 
V(k). To stress the similarity to acoustic phonons, the 
last equation can also be rewritten as 



H = s ^2(gk+TtPkP-k + 9kQkQ-k) 



(B3) 



where 



P k =(4 + a-fc)/\/2 
Qk =i{a\ - 0-fe)/V2 
9k =\V(0) - V(k) + V(0) - V(k)} 



(B4) 



so that canonical commutation relations hold, [Pk, Q q ] = 
iS(k + q). For k ~ we have g k oc fc 2 , g k+v ~ 2V"(0) 
the spectrum w(fc) = 2s^Jg k g k+7r ~ fc is linear for small 
fc. The analogy with acoustic phonons is incomplete, 
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because in the second magnetic Brillouin zone (close to For the nearest neighbor interaction J this formula 
the point k ~ tt) the role of Pk and Qk is reversed with simplifies to: 
respect to the form of their correlation functions. 

For the dynamical susceptibility, x xx i the representa- 
tion in terms of a is the same for both sublattices and we 
have Sfc = y/sPk. The equations of motion read 

^= 2 »«* (B5) S ^TWW- ,B7) 



and hence 



d t Qk = ~ 2sg k+7r Pk 
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